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Abstract 

Hadronic interaction models at cosmic ray (CR) energies are inherently uncertain due to the lack of a 
fundamental theoretical description of soft hadronic and nuclear interactions and the large extrapolation 
required from collider energies to the range of the most energetic cosmic rays observed (> 10^° eV). Model 
uncertainties are evaluated within the QGSJET-II model, by varying some of the crucial parameters in 
the limits allowed by collider data, and between QGSJET-II and other models commonly used in air 
shower simulations. The crucial parameters relate to hard processes, string fragmentation, diffraction 
and baryon production. Results on inelastic cross sections, on secondary particle production and on the 
properties of air showers measured by ground detectors from energies of 10^^ to 10^^ eV are presented. 



1. Introduction 

Direct measurements of the primary cosmic rays (CR) with energies > 10^^ eV are generally difficult 
due to their exceedingly low flux. Instead, their properties are reconstructed from the shape and particle 
content of the extensive air showers (EAS) they produce in the atmosphere. The reconstruction is 
based on numerical models of the air shower development. As there is currently no reliable fundamental 
theory describing soft hadronic interactions at high energies, large systematic uncertainties limit the 
precision of the results. Most current hadronic interaction models use Gribov-Regge theory [l[ of multi- 
Pomeron exchange between nucleons as the basis for the treatment of high-energy, soft interactions, 
which are prevalent in air showers. Perturbative quantum chromo dynamics (pQCD) can describe hard 
interactions with high pj^, which are rare in cosmic ray interactions, but become more important at 
higher energies. In addition, diffractive interactions, collisions of nuclei, and interactions and decays 
of all possible secondary particles at energies from MeV to beyond 10^° eV are needed for a complete 
simulation of an air shower. The generalisation from nucleon-nucleon to hadron-nucleus and nucleus- 
nucleus collisions is usually performed via the Glauber-Gribov formalism 0, taking into account 
inelastic screening and low mass diffraction effects. Particle tracking and electromagnetic interactions 
are straight forward to simulate. The major problem is the seamless and coherent combination of the 
different parts of hadronic models. Different numerical codes implement the same theoretical ideas in 
different ways, with approximate agreement at lower energies where collider data are used for tuning, but 
diverging in the region where extrapolations are required, e.g. to ultra high energies, or to very forward 
emission directions. 

An important aspect of the analysis of experimental data, is the estimation of systematic uncertainties 
of air shower observables, due to imperfections of the current interaction models. The sensitivity of the 
simulated EAS observables to various extrapolations of characteristics of hadronic interactions, such as 
the inelastic cross section, the multiplicity of secondary particles or the relative energy loss of the leading 
(most energetic) hadrons, towards very high energies has been recently investigated in While that 
study is certainly very illustrative, the variations considered can hardly be achieved in realistic models, 
without being in contradiction with available accelerator data or with fundamental principles of the 
quantum field theory. For example, one of the options investigated in Ref. Q assumed the inelastic 
hadron-nucleus cross section to rise asymptotically as In^ s (s being the square of the cm. energy of the 
collision), thus violating the unitarity bound. 

In the present work a different strategy is chosen. We investigate the potential spread in the high 
energy extrapolations of hadronic interaction characteristics and the related uncertainties of air shower 
simulations within a particular model framework. We quantify the uncertainties of the model, by varying 
some crucial model parameters within the limits allowed by accelerator data. In this way, we also take into 
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account the correlations between different interaction characteristics. For example, constraints on the rise 
of the total proton-proton cross section with energy come also from the data on the pseudorapidity (rj) 
density of secondary charged hadrons. The main question addressed here are how uncertain predictions 
can be (within a particular model) and if the uncertainties can cover the differences between currently 
available Monte Carlo generators of hadronic interactions. 

For our study we chose the latest version of the QGSJET model, QGSJET-II-3 ([H, released in 2006), 
because QGSJET describes quite well a large set of experimental air shower data at energies from 10^^ 
to 10^° eV. Additionally, QGSJET-II provides a microscopic treatment of nonlinear interaction effects in 
hadronic and nuclear collisions, described by Pomeron-Pomeron interaction diagrams, which become very 
important at very high energies. We construct a number of alternative parameter sets (options 2 to 6) 
for the model and compare the results obtained with the standard QGSJET-II (option 1) and with two 
models often used in air shower physics, namely SIBYLL 2.1 Q and EPOS 1.99 0. While the original 
SIBYLL was based on the "minijet" idea, ascribing the energy evolution of the interaction cross section 
and particle production solely to the contribution of hard processes (minijet emission), the actual model 
version is rather similar to the Gribov-Regge-type models. A crucial difference between QGSJET-II and 
SIBYLL 2.1 is that parton saturation effects in SIBYLL are modelled by means of an energy-dependent 
cutoff for the minijet production. EPOS employs a phenomenological approach for nonlinear interaction 
effects, while addressing an important problem of energy and momentum correlations between multiple 
scattering processes at the amplitude level EPOS has a significantly larger number of adjustable 
parameters than QGSJET-II, but a substantially larger set of accelerator data has been used for its 
calibration. The distinctive feature of EPOS is the enhanced production of baryons, which leads to more 
energy in the muonic component of showers Q . The standard versions of all these models are available 
in the framework of the latest CORSIKA air shower simulation package ^IQ] . 



2. Parameter Variations in QGSJET-II 



Among the crucial parameters responsible for the high energy extrapolation of the model are those 
related to inelastic diffraction, the relative contributions of soft and hard processes, and the distribution of 
energy and momentum between secondary particles. To investigate the corresponding allowed parameter 
range, five versions (options 2-6) of the QGSJET-II model have been created, with varying parameter 
settings, which are tuned to reproduce (within experimental accuracy) the same set of accelerator data 
as used for the original model calibration. The new data obtained at the Large Hadron Collider (LHC) 
have not yet been used in the analysis. The effects of the parameter variations are then investigated for 
both single interactions and air showers as a whole. Option 1 corresponds to the standard parameter 
settings of QGSJET-II (see Tab. [J). 

Diffractive interactions are very important for the shower development, as they transport the primary 
energy efficiently through the atmosphere. Both, projectile and target nucleons can undergo low-mass 
diffraction excitation. The latter is treated in QGSJET-II in the framework of the Good- Walker-like 
multi-channel eikonal approach 
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|12| |: A hadron is represented by a superposition of a number of 
elastic scattering eigenstates which undergo different absorption, depending on their relative strength Ai 
for the interaction. For the 2-component scheme realised in QGSJET-II, the ratio A1/A2 defines thus the 
strength of diffraction, which rises for more asymmetric superposition (larger A1/A2) and disappears for 
Ai = A2. In addition, a stronger diffraction implies also stronger inelastic screening effects which reduce 
the predicted total and inelastic cross sections. The parameters are varied in two different ways to increase 
the proportion of diffractive events and influence the rise with energy of total inelastic hadron-hadron and 
hadron-nucleus cross sections. In the default version of QGSJET-II (option 1) the corresponding ratios 
are Y'ip) =4 for protons and = 4.7 and 5.7 for pions and kaons respectively. The corresponding low 
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Table 1: Parameter settings of six options of QGSJET-II. Option 1 represents the standard settings of 
QGSJET-II. See text for further explanations. 
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Figure 1: Top: Inelastic cross sections for proton- Nitrogen and 7r-Nitrogen collisions. 
Bottom: Inelasticity of p-Nitrogen and 7r-Nitrogen collisions. 



mass diffraction cross section for the incident hadron changes between 0.9 and 1.6 mb for proton-proton 
and between 1.2 and 1.7 mb for pion-proton interactions in the laboratory energy range of 10^ to 10* 
GeV. In option 2 we enhance low-mass diffraction by choosing ^{p) — 7 and j^(7r) = 9, j^iK) = 12, 
while in option 6 we leave proton diffraction unchanged but decrease diffraction for pions and kaons, 
j^(7r, K) = 3. As an illustration of the discussed modifications, in option 2 the low mass diffraction cross 
section of the projectile hadron at 10^ and 10* GeV is increased to 1.6 and 3.1 mb, respectively, in pp 
collisions and to 2.0 and 3.2 mb in np interactions. 

In the Gribov-Regge framework, QGSJET takes into account contributions of both soft and semihard 
parton cascades. In the latter case, the virtuality (g^) cutoff Qq defines the transition from the non- 
perturbative soft < Qq) to the pcrturbative hard > Qq) part of the cascade. Additionally, in 
QGSJET-II the shadowing and saturation effects are accounted for by the soft (Iqp < Qq) partons. Thus, 
the performance of the model depends strongly on this cutoff parameter and on the parton distribution 
functions obtained from deep inelastic scattering. As the semi-hard contribution rapidly increases with 
energy, one can significantly reduce the rise of the cross sections and the secondary particle production 
with energy when this cut-off is increased from its default value 2.5 GeV'^ to 4 GeV'^, as done in option 
3. 

Options 4 and 5 vary the energy-momentum partition between elementary production processes and 
string fragmentation. In option 4 the "baryon junction" mechanism (BJM) jl^l, related to di-quark 
valence quark interactions, is switched off, leading to a slower energy rise of the inelasticity in proton- 
proton and especially in proton-nucleus collisions. In option 5, in addition the parameter Oq for the 
momentum distribution (~ of constituent (anti-)quarks - string ends (SE) was modified from the 

default value 0.5 to 0.7, which also enhances the leading particle effect in high energy interactions. 

In addition, with the exception of option 4, which is identical to the default model version apart 
from the above-discussed modification, for each new parameter set created a re-tuning of other model 
parameters was performed in order to keep the agreement with accelerator data. The corresponding 
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Figure 2: Left: Average multiplicity of charged particles and 

right: average pj^ for proton- Nitrogen interactions as function of lab energy. 
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Figure 3: Top: dN/drj distributions for charged particles from p-Nitrogen (left) and Nitrogen-Nitrogen 
(right) collisions at E\ah — 10^^ eV. In these plots lines from options 1, 4 and 6 lie on top of one 
another. Bottom: dE/drj distributions (left) and Feynman-x distributions (right) for charged particles 
from p-Nitrogen collisions. All plots are for 10^^ eV. Colours as for Fig. [51 



changes were minimal for option 6 (weaker diffraction for pions and kaons), which concerned the Pomeron- 
pion interaction vertex only. For option 5, a significant modification of the string fragmentation procedure 
was necessary to compensate the decrease of particle multiplicity due to shorter strings. An opposite 
procedure has been used for option 2, where the particle yield from fragmenting strings was decreased. 
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For options 2, 3, and 5, most of the basic model parameters, like Pomeron intercept and slope, Pomeron- 
hadron couplings, and parton (quark and gluon) distributions in the Pomeron had to be re-tuned, which 
was necessary to remain in agreement with the data on total and elastic hadron-proton cross sections, 
elastic scattering slopes, and on the proton structure function F2. 

The main parameters modified in the different model versions are listed in Tab. [1] While the 
parameters investigated here are certainly crucial ones, it is not clear whether they are the only ones 
that matter. The comparison with EPOS shows that there are other mechanisms which can influence 
significantly the shower properties and consequently affect the interpretation of air shower experiments. 

3. Single Interactions 

In this section the particle production in individual collisions is investigated. It relates directly to the 
particle content and shape of showers, and is crucial for energy determination and composition analysis 
in air shower experiments. Overall, for all options a reasonably good agreement between the available 
collider data at low energies and the simulations was achieved. The evolution at high energies Fig. [1] 
shows the inelastic cross sections and the inelasticity of p- Nitrogen and 7r-Nitrogen collisions as a function 
of energy, for all the models. The cross sections of the QGSJET variants have similar shapes and the 
differences between the options are within 10%, even at the highest energies, with the smallest values 
being obtained for options 2 and 3 (higher diffraction and higher Qo-cutoff, respectively). EPOS shows a 
cross section very similar to the QGSJET models, which could be expected as both models are genuinely 
based on Gribov-Regge theory. However, SIBYLL shows a rounder shape with a much steeper increase at 
high energies, most likely being related to the adopted parametrisation for the energy-dependent cutoff 
of minijet production. 

The inelasticity is the fraction of energy used for the production of secondary particles, i.e. 1 — 
Eip/Etot, where Eip is the energy of the leading baryon for proton collisions and the leading charged 
meson for pion collisions. The evolution of the inelasticity as function of -Eiab is shown at the bottom 
of Fig. [H These plots show that all QGSJET model variants produce a similarly shaped dependence of 
inelasticity with energy, with a characteristic minimum around 10"^ — 10^ GeV, but some are shifted to 
lower inelasticity values by up to 0.05. This is also true for the inelasticity of pion-Nitrogen collisions, 
except that the shifts are much smaller. Both, EPOS and SIBYLL, show significantly different shapes in 
inelasticity vs. energy, but with similar values as QGSJET at around 10^^ eV and 10^^ eV. 

Fig. [2] shows the average number of charged particles produced in p-Nitrogen collisions, {rich) is very 
similar for all the models, up to a lab energy of about 10^^ eV, where the models start to diverge. At 
the highest energies the QGSJET variants all produce 2-5 x more secondaries than SIBYLL or EPOS. 
All models predict a similar average pj^ growing with energy, only EPOS gives smaller p_L values below 
10"* GeV and larger ones above 10* GeV (see Fig. [2l right). 

The large differences in particle numbers can also be seen in the pseudorapidity distributions of 
charged particles (see Fig. |3]). QGSJET produces the largest numbers of secondaries, and there are big 
differences among its variants. At 10^^ eV QGSJET option 5 produces almost 50% more, and option 3 
only half the number of particles than option 1 in the central region (77 « 0). SIBYLL and EPOS are even 
lower than option 3. Recent measurements at LHC in the central pseudorapidity range indicate indeed 
that QGS JET-II tends to overestimate the central pseudorapidity density of secondary hadrons (see 
for a detailed comparison. However, most of these particles are emitted with low energies, compared 
to the energy of the parent hadron. Although differences can be also seen at larger pseudorapidities, 
these are of much smaller magnitude than in the centre. While being relevant for calculations of the 
niuon content of air showers, the hadrons in the central rapidity region are of little importance for the 
longitudinal shower development. The fact that the forward particles are more important can be seen 
from the distribution of the energy flux (dEi/drj) vs. pseudorapidity in Fig. |3]bottom, left. In the forward 
region, where diffractive interactions contribute and where most of the interaction energy is going, all the 
models are fairly close together. Very little energy emitted in the central region, despite the large number 
of particles. EPOS shows good agreement with the QGSJet variants. SIBYLL, however, has prominent 
peaks in dE/drj in the forward and backward regions. Also the Feynman-x distributions show that there 
is fair agreement between all the models in the very forward region (x-p — PL/Pmax ^ 0.95). 

The 77 distributions also allow assessment of the treatment of nuclear interactions in the models, 
which is important for nuclear composition studies of cosmic rays. As Nitrogen-Nitrogen collisions are 
perfectly symmetric, also the j] distribution should be symmetric. In QGSJET nuclear interactions are 
modelled as two clouds of nucleons penetrating each other and interacting. As the modeling is done in a 
symmetric way also the outcome is symmetric. Also EPOS produces a symmetric distribution. SIBYLL 
models nuclear interactions with the semi-superposition model which treats a nuclear projectile as a 
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superposition of individual nucleons and a nucleus-nucleus collision as a number of independent nucleon- 
nucleus collisions. This is by construction asymmetric and therefore the Nitrogen-Nitrogen collisions 
result in asymmetric pseudorapidity distributions. SIBYLL produces > 4x more particles in Nitrogen- 
Nitrogen than in p-Nitrogen collisions, while for QGSJET and EPOS this factor is only about 2. Thus, 
for composition studies based on differences in the muon yield in showers from different nuclei, it is better 
to rely on QGSJET or EPOS. 

The large differences in ry between models and the comparison with first collider data at high energies 
indicate that some aspects of the models are not correct and need adaptation. Although currently there 
are only few measurements in the relevant forward region and LHC is not yet running at its maximum 
energy, LHC results, in particular, the results of the LHCf experiment [l^ provide already constraints 
to some of the model parameters. Future measurements of the total proton-proton cross section and of 
particle production at -y/s = 14 TeV will greatly help to tune the models further. 

4. Air Showers 

Showers were simulated at energies of 10^^ , 10^^ and 10^^ eV, i.e. at typical energies relevant to 
Cherenkov telescopes (such as HES S VERITAS [17} and MAGIC [ll]), small air shower arrays (such 
as KASCADE [19| and LHAASO [20(|), and experiments for the highest energies (e.g. the Telescope 
Array [2l| or the Pierre Auger Observatory f22'|), respectively. The first interaction point of the primary 
particle was fixed at typical values for the respective primary energy, to exclude the effect of shower- 
to-shower fluctuations due to the varying point of first interaction and to concentrate on differences in 
average shower shape. For each energy the average longitudinal development and the lateral distribution 
of particles at ground level have been determined. Also characteristic observables were evaluated, such as 
the average values of the atmospheric depth of the shower maximum, X^^^, the particle number at the 
shower maximum, A^max, the Cherenkov photon lateral distribution at ground level, the electron-to-muon 
ratio at ground level, the energy that would be released in an Auger- like water Cherenkov detector at 
1000 m core distances (S'looo) |23| and the time in which the integrated signal grows from 10% to 50%, 
ti/2 [23]. Simulations have been performed for the QGSJET variants as well as for SIBYLL 2.1 and EPOS 
1.99 to allow comparisons between the scale of the uncertainties due to parameter variations within a 
model, and due to model-to-model variations. 

4-1. Low Energies (10^^ eV) 

At 10^^ eV very little variation is found in both the lateral distribution and the longitudinal develop- 
ment, for the QGSJET options, with differences being < 10% for all quantities investigated. This is not 
surprising as the energy is close to the accelerator energies at which the models have been tuned. 

However, as the number of shower particles reaching ground level at this energy is very low, a more 
appropriate observable is the lateral distribution of Cherenkov photons at ground level. The number 




Figure 4: Left: Percentage difference of the Cherenkov photon density p (for A = 300-600 nm) compared 
to the standard QGSJET-II model (i.e. Pi/pstd) — 1) at 2000 m altitude. Results from 10* vertical proton 
showers of 10^^ eV. Atmospheric absorption is accounted for. Error bars show statistical errors of the 
mean values. Right: Fraction of total interaction energy emitted in photons as a function of energy 
(Ej/Etot) 
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Table 2: Simulated mean observables for proton showers of energy 10^^ eV at 0°. The errors given are 
the statistical errors of the mean values. Note that the position of the first interaction was fixed at a 
depth of 33 g/cm^. 



of Cherenkov photons relates to the primary energy, but is different for different primaries. Proton 
induced air showers are the main background in gamma ray experiments using the atmospheric Cherenkov 
technique and their simulated appearance depends on the hadronic interaction models. Small differences 
are seen in the lateral distributions (Fig. 2]), resulting from the differing fractions of energy that is 
channeled into the electromagnetic (EM) component. This is visible in QGSJET option 2, where the 
biggest fraction of energy is going into the EM component, leading to the largest number of Cherenkov 
photons and lowest number of muons. SIBYLL and EPOS predict the largest numbers of Cherenkov 
photons, with about 7% more than the standard QGSJET option. This difference alone contributes 7% 
to the systematic energy error of proton showers in Cherenkov telescope experiments. The systematic 
shift is lower than the energy resolution of current instruments (w 15% for 1 TeV photons), so it should 
not affect much measurements of the background in ACT instruments. 

The ability to identify electromagnetic gamma ray showers from the overwhelming background of 
hadronic events directly determines the sensitivity of a Cherenkov telescope. It is based on the analysis 
of the shape of the shower image and depends on hadronic models, too. Hadronic showers produce 
more diffuse and larger images in the camera of a Cherenkov telescope than compact electromagnetic 
showers. Those hadron showers that are misidentified as gamma showers are usuall y ch aracterised by a 



large fraction of the shower energy being dumped early on into the EM component 25 1. 

The EM fraction in Fig. |4] looks similar to the inelasticity shown in Fig. [1] with values of about 1/3 
of the ones of the inelasticity, as one would expect if 1/3 of the pious produced are 7r°s. Again, most 
QGSJET variants produce a lower EM fraction than the standard QGSJET model, and SIBYLL and 
EPOS show different shapes. Overall, SIBYLL and EPOS were found to produce a larger EM fraction 
than QGSJET above 10"'^^ eV. A larger EM fraction in the first interaction means a more photon-like 
appearance of the shower and a higher chance to misidentify a proton-induced shower as a gamma 
ray. Thus, all QGSJET variants give an equal or smaller rate of photon-like events than the QGSJET 
standard. Both, SIBYLL and EPOS will produce more gamma-like proton showers, but to quantify this 
would require full simulation of the telescope array and the reconstruction procedure. 

This difference in the EM fraction can be partly explained by the differences in the inelasticity of the 
models, changing the amount of energy available for particle production. Once this effect was accounted 
for the QGSJET variants were found to give here a constant value of about 31% whereas SIBYLL is 
about 1% higher and EPOS 1% lower. As EPOS leaves more energy in secondary baryons, it is not 
surprising that the fraction in tt'^s, and thus photons, is smaller. 

4-2. Medium Energies (10^^ eV) 

At medium energies (^ 10^^ eV) detector arrays usually measure shower particles arriving at ground 
level, with some arrays being able to distinguish between e^, photons and muons. Therefore, Tab. [3] 
gives some typical shower observables for the different models. The left side of Fig. [5] shows differences 
in electron and muon numbers as function of core distance (lateral distributions) and ionisation energy 
deposit as a function of atmospheric depth (longitudinal distributions) for 10^^ eV proton showers. At 
10^^ eV differences between the QGSJET model variations are less than ±10% for electrons and muons 
at core distances of 10 to 1000 m (the most relevant core distances at these energies are < 200 m). 
The longitudinal distributions show very good agreement around the shower maximum (~550 g/cm^) 
and differences up to ±10% at ground level (here, 880 g/cm^). SIBYLL lies well within the spread 
of the QGSJET options, but EPOS shows much flatter lateral distributions, for both electrons and 
muons, with many more particles at large distances from the shower core and overall about 20% more 
muons. At ground level EPOS has 15% fewer electrons and 15% more muons than the reference model. 
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Figure 5: Top: Average percentage difference (compared to standard QGSJET-ff) of the particle number 
N as function of the core distance at 1452 m ahitude (i.e. A^^/A'std — 1) for electrons and muons. Bottom: 
Average percentage difference (compared to standard QGSJET-II) of the ionisation energy deposit as a 
function of the atmospheric depth for the em component and for muons. Left: 300 vertical proton 
showers of 10^^ eV. Right: 100 proton showers at 20° with 10^^ eV. Shaded error bands show statistical 
errors of the mean values. 
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630 


±2 



Table 3: Simulated mean observables for proton showers of energy 10"'^^ eV at 20°. The errors are the 
statistical errors of the mean values. Note that the position of the first interaction was fixed at a depth 
of 38 g/cm^. 



Consequently, all energy assignments would be different on the 15-20% level and composition analysis 
were greatly affected by these differences. They are big enough for the KASCADE experiment to tell the 
difference. Indeed, there is evidence from KASCADE data |26| that the primary proton spectrum agrees 
better with direct measurements when QGSJET and SIBYLL are used for the analysis, than when EPOS 
is used. 

The lateral distribution of QGSJET option 2 (enhanced diffraction) is very similar to that of SIBYLL 
at both 10"^^ and 10"^^ eV, suggesting that one of the major differences between the two models is the 
rate of diffraction. However differences can be seen between option 2 and SIBYLL in the results of single 
particle interactions, indicating that there are other important physical distinctions. 



4.3. Ultra High Energies (10^^ eV) 

At 10^^ eV the detector arrays must be sparse and no longer the total particle number is measured, 
but the particle density as a function of distance from the core. The energy is then determined from the 
density at some distance 27|. The time structure of the shower front, measured in the array detectors, 
carries information on the shower development and the primary mass. With fluorescence telescopes the 
longitudinal profile and the position of the shower maximum can be directly measured, from which the 
energy and the composition can be inferred. Commonly used observables at 10^^ eV are listed in Tab. [31 

As expected, the deviations from standard QGSJET are largest at this energy, as the extrapolations 
of model parameters from collider energies are largest. Deviations in the lateral distributions of electrons 
and photons are less than +10% within 1000 m from the shower axis. For muons the variations within 
the QGSJET models are less than 20% at 1 km core distance, but all the options produce always smaller 
muon numbers than the standard QGSJET. SIBYLL gives almost the same results for electron and muon 
lateral distributions as QGSJET option 2, the version with the enhanced proton diffraction. Again, EPOS 
is distinctly different, with 25% more muons at ground level for all core distances. The QGSJET and 
SIBYLL muon numbers at ground are within to -15% of the standard version. The EPOS electron 
numbers at ground are similar to those of the other models, owing to the fact that for the angle chosen 
the shower maximum is close to the ground. For more inclined showers larger differences would occur. 

In the longitudinal distributions, the QGSJET options 2, 3 and 4 and SIBYLL have markedly lower 
particle numbers and energy deposits (40-20%) in the upper part of the shower than the standard. EPOS 
shows a similar trend, again developing slowly early in the shower, but quickly producing larger numbers 
of muons, with almost 30% more muons at ground level. 

To reproduce data from experiments directly measuring muon numbers such as KASCADE and HiRes- 
MIA and from the Pierre Auger Observatory, where the muon content can be inferred indirectly |28| . 
more muons in simulations seem to be required. EPOS seems to address the muon deficit problem, but 
would introduce a significantly different energy scale. 

Significant differences are also seen in S'loooj a common observable used for the energy reconstruction 
in the 10^^ eV region (see Tab. [S]). The energy reconstruction based on option 2 would give ^-^10% 
higher cosmic ray energies than the one using the standard QGSJET option, while EPOS would give 16% 
smaller energies. The differences in average ATmax for the QGSJET models are less than 30 g/cm^ (4%), 
with all the models producing deeper Amax values than standard QGSJET. The differences in ti/2, which 
is used for composition analysis, are less than 5%, with SIBYLL and EPOS giving the extreme values. 
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5. Conclusion 



The analysis of sliower simulations demonstrates the increase in systematic uncertainties (introduced 
by variation of model parameters) as the primary particle energy increases, and with it the extrapolation 
of the models. Our study indicates that parameter variations within one model generally cannot reproduce 
the full variance possible between models. None of the parameter variations investigated in QGSJET was 
able to mimic the behaviour of the EPOS model: At any energy, EPOS has far more muons than all other 
models and specifically than the QGSJET standard version. Differences in several observables amount to 
< 10% in the TeV range and to 20-30% in the PeV and EeV ranges, causing sizeable systematic errors on 
the energy scale and composition analysis. Within QGSJET, at all primary energies the largest effects 
are observed due to the increase of the rate of diffraction for protons and pions, making diffraction the 
feature with the largest impact on the overall development of the shower. 

EPOS is a relatively new model and still has to prove its qualities by being compared in detail to 
results from different (past and present) experiments from TeV to EeV energies. 

The lateral and longitudinal distribution pattern of the electromagnetic and muonic components in 
air showers shows complex variations between models, which depend on the energy and the level of 
shower development, the zenith angle, etc. Therefore, complete simulations of showers, detectors and 
analysis procedures, using different interaction models are needed for a meaningful data analysis. The 
model uncertainties at high energies are of a size that makes analysis of nuclear composition at >EeV 
energies rather difficult (with the observables considered here), as differences between proton and iron 
induced showers are not much larger than the uncertainties due to the models. On the other hand, the 
uncertainties are large enough that current experiments are able to tell which model fits best. Also, new 
data from LHC and RHIC do increasingly constrain the models, and help to make the extrapolation 
to highest energies more reliable. But ultimately, only a coherent description of cosmic ray phenomena 
and hadronic interaction physics over the full energy range will be a convincing proof of the correct 
interpretation of cosmic ray data. 
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